Precession of Vortices in Dilute Bose-Einstein Condensates at Finite Temperature 
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We demonstrate that the precessional frequencies of vortices in Bose Einstein condensates (BECs) 
are determined by a conservation law, and not by the lowest lying excitation energy mode. We 
determine the precessional frequency for a single off-axis vortex and vortex lattices in BECs using 
the continuity equation, and solve this self-consistently with the time-independent Hartree-Fock- 
Bogoliubov (HFB) equations in the rotating frame. We find agreement with zero temperature 
calculations (Bogoliubov approximation), and a smooth variation in the precession frequency as the 
temperature is increased. Time-dependent solutions confirm the validity of these predictions. 
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Quantised vortices are perhaps the single most striking 
manifestation of superfluidity. The dynamics of vortices 
in Bose Einstein condensates (BECs) have been of par- 
ticular interest [1] with earlier theoretical work [2, 3,0,[1| 
concentrating on zero temperature, neglecting the ef- 
fects of the non-condensate on the overall dynamics of 
the vortex. The Bogoliubov spectrum in this case has a 
lowest-lying energy that is negative, the anomalous mode 
B El El- In the frame rotating at the frequency corre- 
sponding to this energy, the lowest-lying energy in the 
Bogoliubov excitation spectrum is zero, leading to the 
conclusion that this mode frequency corresponds to the 
precessional frequency of the vortex 0, 0, [B[ . 

To go to finite temperature, one can introduce the 
Hartree-Fock-Bogoliubov (HFB) treatment @. Solving 
the HFB equations in the laboratory frame for an on- 
axis vortex at finite temperature in the so-called Popov 
approximation [6] yields a positive lowest lying energy 
mode, referred to as the lowest core localised state 
(LCLS)[7|, M. This treatment was generalised to an off- 
axis vortex [9[ and provided a lower bound for the LCLS 
energy. In both cases it is argued that the thermal cloud 
acts as an effective potential, thereby stabilising the vor- 
tex, but both fail to take into account the dynamics of the 
thermal cloud itself. The association of the LCLS energy 
with the precessional frequency of the vortex Q led to the 
conclusion that the vortex precesses in the direction op- 
posite to the condensate flow around the core. This is in- 
consistent with experiment [10[ and also seems intuitively 
unreasonable, since zero temperature predictions would 
suggest otherwise d, 0, H| and one might reasonably ex- 
pect continuous variation with temperature. However, 
in later work on off- axis vortices [11], Isoshima et. al. 
conclude that the precessional frequencies and LCLS en- 
ergies are indeed uncorrelated, but assert rather that the 
LCLS energy is responsible for the inward/outward mo- 
tion of vortices. They have no explicit equation for the 
prediction of precessional frequencies, but rely on the as- 
sertion that the LCLS energy depends linearly on the 
rotational frequency, and solve the problem by iteration. 

Here we make use of the continuity equation to make 
an a priori prediction of the vortex precessional fre- 
quency and solve the HFB equations self-consistently in 



the frame rotating at this frequency. This equation is also 
valid for vortex lattices, and stationary solutions can be 
found in the frame rotating at this precessional frequency, 
provided the vortex interactions are not too strong. 

In a frame rotating with angular frequency ^, the time- 
dependent HFB formalism [12] consists of the generalised 
Gross-Pitaevskii equation (GGPE) for our condensate 
wavefunction 

ih£$(r,t) = (hn(r) + (7(n c (r,t) + 2n(r,t)))$(r,t), 
+<7m(r,f)**(r,t) 

(1) 

and the Bogoliubov-de Gennes equations (BdGEs) 



with 
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Mr) = -£v 2 +m(rxV) + y T (r), 

£(r, t) = h(r) + 2g (n c (r, t) + n(r, t)) - /i, (3) 

M(r,t) = #($>,£) +m(r,£)). 

Here g = 47rh 2 a s /m is the interaction strength, a s the s- 
wave scattering length, m the atomic mass. The conden- 
sate density n c , thermal density ft and anomalous density 
rh have their usual definition Q. 

Defining the current density j(r, t) = 
-if ($*(r,t)V$(r,t) -$(r,t)V$*(r,t)), one obtains 



from (pQ) the continuity equation 



\$(r,t)\ 2 + £Vj(r,t) = m.(r x V) |$M)| 2 
-ig (m(r,i)$* 2 (r,t) - m*(r, t)$ 2 (r, i)) . 



(4) 



The corresponding time-independent HFB equations are 
given by [6] 

M$(r) = (h n (r) +g(n c (r) + 2n(r))) $(r) + <?m(r)$*(r), 

(5) 
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where /i is the chemical potential, and e q are the quasi- 
particle energies. 

We note, however, that the quasi-particle amplitudes 
calculated using this formalism are not, in general, or- 
thogonal to the condensate [13[. This can cause anoma- 
lous predictions for Q in various temperature regimes. 
We can derive a suitable orthogonal formalism, i.e. one 
where J dr$*(r,t)u q (r,t) = f dr$(r,t)v q (r,t) = as 
follows: We write the field operator [14] as ^(r,i) = 
a c (t)(j)(r,i) +r)(r,£), where a c (t) = f dn/>(r, t)0*(r, t) 
and a\(t) are the annihilation/creation operators cor- 
responding to the condensate. 0(r, t) is the conden- 
sate wavefunction normalised to unity such that = 
y/Nc(f> where N c is the number of condensate atoms, and 
f)(r,i) = j dr'Q(r, r')i/;(r', t) is the fluctuation opera- 
tor corresponding to the thermal cloud, with Q(r, r') = 
<5(r, r') — 0(r, t)(j)*(r f , t). Using the commutation rela- 
tions [a c (t),al(t)] = 1, [r)(r,t),7)t( r ^)] = Q(r,r') and 
[a c (t), tf(r, t)~\ = 0, and applying the Bogoliubov trans- 
formation and the usual mean-field approximations [6], 
one obtains the orthogonal time-independent HFB equa- 
tions consisting of the GGPE ([5]) and the orthogonal 
BdGEs given by (|6j), but where the Cq(t) and Ai(r) 
operators are replaced by 



Ai(r,r') - Jdr'Q(r,r')£ Q (r') 
M{ry) -> fdr'Q(r,r')M(r') 



(7) 



This ensures the orthogonality of the condensate and 
thermal states. In view of this orthogonality, making 
the assignment uq — > </>, —> — </>* satisfies the BdGEs 
([6]) with a zero energy. Thus the lowest excitation en- 
ergy is zero (which is certainly not the case for ordi- 
nary HFB), with quasi-particle amplitudes proportional 
to the ground state, thereby satisfying Goldstone's theo- 
rem. However the effect on higher energy modes is neg- 
ligible, and the energy gap [15] still remains an issue. 

All time-independent calculations presented here were 
performed using this orthogonal HFB. The results with- 
out the orthogonalisation are essentially the same, ap- 
part from some slightly anomalous behaviour in Q for 
vortices close to the axis at low temperature, where the 
orthogonality of the solutions becomes important. We 
note that the Popov approximation, as opposed to HFB, 
violates particle conservation, and is therefore not a suit- 
able formalism for time-dependent simulations. HFB can 
also be derived from a variational standpoint, and rep- 
resent a "conserving" approximation [6[. It is from this 
basis that we use HFB for our simulations. As a valid- 
ity check, time dependent simulations where the BEC is 
initially displaced from the axis yield Kohn mode [l6[ 
oscillations at precisely the radial trapping frequency. 

To proceed, we write the continuity equation in terms 
of the real and imaginary parts of the condensate wave- 
function. Defining R(r) = Re($(r)) and J(r) = 
Im($(r)), and noting that in the time-independent case 



§- t |$(r,£)| 2 = 0, we obtain 
n. (R(r x V)/ + I(r x V)R) = ^ (RV 2 I - IV 2 R) 

(8) 

Our model consists of a dilute Bose gas consisting of 
2000 87 Rb atoms in an axially symmetric harmonic trap, 
with radial and axial harmonic trapping frequencies of 
Vt r = 2tt x 10Hz and fl z = 2tt x 400Hz respectively. Thus 
the axial confinement is sufficiently strong that all excited 
axial states may be neglected. Hence the BEC may be 
treated as two-dimensional from a computational view- 
point, although the trap geometry is not critical for our 
conclusions. Thus, restricting ourselves to two dimen- 
sions, and defining 
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and integrating over all space, we find that 

n= / A(r,9)B(r,0)rdrd0/ / / r 2 (A(r, 0)f rdrdO. 
Jo Jo Jo Jo 

(9) 

We can now solve equations ([5]), (j6]) and ([9]) self- 
consistently in the frame rotating at angular frequency 
Q to find stationary solutions for a precessing vortex. A 
vortex at position (ri,#i) may be created by expanding 
the condensate wave function in terms of modified La- 
guerre basis functions |x^( r '^)} as f°ll° ws: 



(i) 



*(r,fl) = ^;a in x}; ) (r,fl), 



(10) 



provided (n, 6\) is not a root of & n (r, 0), where we define 
the computational basis functions 



X^(r,f)=6n(r,< 
and where 



6n(n,6>i) 
£io(nA) 



CioM), (11) 



no 



2n\ 



V2tt \{n + 



1/2 



- 2 /V^|^l(r 2 ), (12) 



with Ln (x) a modified Laguerre polynomial of order 
n. The superscript (1) in the summation means ex- 
clude (n, I) = (0,1), and indicates a single vortex. 
We derive ([TO]) and (JTTj) by expanding the condensate 
wave function in terms of the complete Laguerre basis 
{£/n( r > n = 0, . . . , I = 0, ±1, . . .} and considering that 

^(n^l) = E^n^n6n(ri^l)=0. 
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FIG. 1: (colour online) (a) Popov condensate (solid line) and 
thermal densities (dashed line) in the plane passing through 
the vortex, and (insert) condensate fraction for the ideal gas 
(dashed line), HFB (circles) and Popov calculations for con- 
densate with Bogoliubov mode cutoff of 100 (triangles), and 
1639 (squares), (b) HFB condensate (solid line) and thermal 
densities (dashed line) in the plane passing through the vor- 
tex, (c) Precessional frequency Q versus vortex position a for 
HFB at temperatures OnK, 2nK, 5nK and 7.5nK (triangles, 
diamonds, circles, squares), (d) Precessional frequency Q ver- 
sus temperature T for HFB at vortex positions 0.1, 0.5, 0.8 
and 1.1 (triangles, diamonds, circles, squares). All units are 
given in radial harmonic oscillator units ro = \fh/mVL r for 
distance and Eq — h£l r /2 for energies and frequencies. 



This procedure may be extended to N v vortices 
#1 ) •> - - - •> ( r N v , Qn v )} by an iterative process, where 
we write 



$(r, 0) 



(Nv) 
In 



(13) 



where we have defined 



{N v - 



X Nv ( r N v , 



(14) 

Stationary solutions are found in the single vortex case, 
and in the multiple vortex case, provided the interactions 
between the vortices are not too strong. 

Figures 1(a) and 1(b) show the condensate and non- 
condensate density profiles in the plane passing through 
the vortex for the Popov and HFB cases respectively for 
a vortex situated at a = 0.5 radial harmonic oscillator 
units from the axis, at temperature T c = 5nK. Note in 
particular that the thermal density for HFB is consider- 
ably less than for Popov. This is due to the upward shift 
in the lower- lying excitation energies. 



One would expect the precessional frequency ft to vary 
smoothly with temperature T, and the time-independent 
off-axis calculations for a variety of different vortex po- 
sitions for the Popov and HFB cases reveal that this is 
indeed the case, also being consist with off-axis GPE cal- 
culations, which coincide almost exactly with the HFB 
T = calculations (see figures 1(c), 1(d)). This is in 
stark contrast with previous @, @, 0] inferences regarding 
vortex precession based upon on-axis HFB-Popov calcu- 
lations with static thermal clouds. These calculations 
therefore resolve this conflict. 

One would also expect the precessional frequency Q 
to vary smoothly with a, so to check the consistency of 
the on-axis and off-axis calculations, one needs to deter- 
mine ft for the on-axis case by extrapolation. Solving 
for the off-axis case in the frame rotating at the extrap- 
olated precessional frequency for the Popov approxima- 
tion, using ©, one finds consistency in the LCLS en- 
ergies. These energies also correspond to the energies 
calculated perturbatively in [ItJ with which we find very 
good agreement at low temperature. One also finds ab- 
solutely no correlation between the precessional frequen- 
cies and the LCLS energies in agreement with [11]. The 
excitation due to clcls is given by the mode density 
corresponding to the LCLS mode, and has nothing to do 
with precession. The precessional frequencies are deter- 
mined by the conservation of flow (the continuity equa- 
tion), and hence by equation (j9]). These calculations also 
predict precession in the same sense as the circulation, 
as one would expect, contrary to [8]. It has been sug- 
gested that the vortex precession inferred from previous 
calculations might be associated with the nutation of the 
vortex about the artificially static effective pinning po- 
tential due to interactions of the condensate with the 
thermal cloud. To investigate this, time dependent simu- 
lations were performed with the thermal density initially 
rotating with the core at the precessional frequency, but 
where the rotational frequency of the thermal cloud is 
gradually reduced. One finds that the vortex core ini- 
tially slows down with the pinning potential to a certain 
point, depending on the strength of the pinning poten- 
tial of thermal density. In principle, given a sufficiently 
strong pinning potential, one should be able to stop the 
precession comletely, though to date this has not been 
achieved. No mutation of the vortex (i.e. precession 
of the vortex about the pinning potential) has been ob- 
served in any of the simulations. We are left to con- 
clude that the LCLS energies play no part in precession 
or nutation of the vortex and the LCLS corresponds to a 
collective response of the gas unassociated with the pre- 
cession of the vortex. We note further that clcls varies 
continuously with T, and we see that €lcls\t^o a l wavs 
takes on some small positive value, associated with the 
quantum depletion. 

We also find consistency in the condensate fractions 
between the on-axis and off- axis calculations, and the in- 
sert in Figure 1 (a) shows a plot of the condensate fraction 
N c versus T(nK) for the off- axis case for HFB, and the 
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FIG. 2: Vortex precession in HFB time- dependent simulations 
showing (a) condensate density, (b) thermal density, (c) con- 
densate phase for a single vortex at a = 0.5, and (d) HFB con- 
densate density, (e) thermal density, (f) condensate phase for 
a triangular vortex lattice for vortices at a — 1.65. All units 
are given in radial harmonic oscillator units ro = yJh/mQ, r 
for distance and E = Ml r /2 for energies and frequencies. 

on-axis case for Popov for a mode cutoff of 100. Ad- 
ditional points are shown for Popov for mode cutoff of 
1639. Shown also is N c versus T for the ideal gas with 
T c « 16.7nK . The results reveal that the mode cut- 
off of 100 is sufficient for T < 7.5nk, but the results for 
T > T c /2 become unreliable. At these temperatures the 
HFB formalism is unreliable in this regime anyway [18[, 
so a cutoff of 100 modes represents an adequate com- 
putational investment. At this cutoff, for T « 5nK, one 
estimates an error of ~ 10% in thermal population which 
does not qualitatively affect our results. 

Time-dependent simulations were performed using full 
HFB for T = 5nK for a vortex at a — 0.5 radial har- 
monic oscillator units from the axis, with a Bogoliubov 
mode cutoff of 100, using the time-independent initial 



state. In all the simulations, the vortex precession de- 
scribes a circle. There is no evidence of dissipation dur- 
ing the time of the simulation of 10 trap periods. The 
precessional frequency Qls — 0.3794^ r was estimated 
using least squares, in good agreement with the value 
of ft = 0.3727£7 r , as predicted in the time-independent 
calculations. The error in the prediction scales as the 
number of computational basis functions which, for prac- 
tical reasons, is only 209 in these calculations. For 839 
computational basis functions, the time-independent cal- 
culations predict a value of Q = 0.3761£l r . 

Figures 2(a),(b),(c) show respectively contour plots for 
the condensate density, the thermal density, and the con- 
densate phase in the x-y plane. We note the phase discon- 
tinuity in figure 2(c) indicating a singly charged vortex. 

Time-dependent HFB simulations were also done for 
a triangular vortex lattice at T = 5nK, with vortices 
symmetrically positioned at a = 1.65 radial harmonic os- 
cillator units from the axis, again using time-independent 
calculations for the initial state. For this symmetric case, 
the vortex precession again describes a circle, and there is 
no evidence of dissipation during the time of the simula- 
tion of 10 trap periods. Figures 2(d),(e),(f) show respec- 
tively contour plots for the condensate density, the ther- 
mal density, and the condensate phase in the x-y plane. 
We note the phase discontinuity in figure 2(f) indicating 
three singly charged vortices. Non-equilibrium studies of 
vortex lattices using the time-dependent HFB formalism 
will form the basis of future investigation. 

In conclusion, we have performed HFB calculations for 
an off-axis vortex using the continuity equation to pre- 
dict the precessional frequencies, which we found to be 
uncorrelated with the LCLS energies, in agreement with 
pj| . The sense and frequency of the vortex precession ex- 
trapolate smoothly from the zero-temperature GPE be- 
haviour removing all previous anomalies. These results 
were confirmed using time-dependent HFB simulations. 

This work was funded through the New Economy 
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Technologies. 
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